{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cea953ad",
   "metadata": {},
   "source": [
    "# 单量子比特门分解算法\n",
    "\n",
    "## 1 课题背景：\n",
    "\n",
    "在容错量子计算中，减少量子过程产生误差的一种常用方法是采用高精度基本门集合来实现任意门操作。典型的基本门集合包括H门、T门、CNOT门。理论研究表明，任何单量子比特门，给定一个特定的精度，总是可以表示为一个有限长度的$H,T$门序列。反之，可以通过为序列设置一个固定的长度来逼近目标量子门，随着序列长度的增加，逼近精度提高。\n",
    "理论上，蛮力搜索总是可以找到任何单量子位门分解的最优解（计算搜索空间中每个可能的门序列与目标门的距离），但其计算复杂度随着序列长度的增加呈指数增长，在实际应用中难以实现。为了应对这一挑战，研究者们已经提出了多种方法，包括SK算法，中间相遇算法，KD树算法等。\n",
    "\n",
    "## 2 项目任务及实现方案：\n",
    "\n",
    "本项目的核心任务是实现一种使用启发式算法高效分解任意的单量子比特门的方法，展现出整体求解时间或者求解返回的基础量子门数量相比之前的经典算法具有不小于一个数量级以上的优势。\n",
    "\n",
    "### SK算法\n",
    "\n",
    "本文代码实现了近似单比特量子门的 Solovay-Kitaev 算法（SK 算法），其目标是将任意单比特 $SU(2)$ 单位门，逼近为基本门集合（$H$ 和 $T$ 门）的组合。\n",
    "\n",
    "1.导入相关依赖"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "28f78841",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "from random import randint, seed\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from utils import helper\n",
    "from utils import ops\n",
    "from utils import state\n",
    "import math\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11df6f4a",
   "metadata": {},
   "source": [
    "ops、helper、state：辅助模块，定义了量子门、状态和辅助函数。\n",
    "\n",
    "2.工具函数定义"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3bd9fdd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_su2(u):\n",
    "\n",
    "  det = u[0][0] * u[1][1] - u[0][1] * u[1][0]\n",
    "  return np.sqrt(1 / det) * u\n",
    "\n",
    "def F_dist(u, v):\n",
    "\n",
    "  return math.sqrt(1 - 0.5 * abs(np.trace(u.adjoint() @ v))\n",
    "\n",
    "def generate_su2(s):\n",
    "  seed(s)\n",
    "  theta =[randint(1,10) for i in range(3)]\n",
    "  t =np.linalg.norm(theta)\n",
    "  theta[0] = theta[0]/t\n",
    "  theta[1] = theta[1]/t\n",
    "  theta[2] = theta[2]/t\n",
    "  su2 =  ops.Rotation(theta,t,'s')\n",
    "  return su2\n",
    "\n",
    "def create_unitaries(base, limit):\n",
    "\n",
    "  gate_list = []\n",
    "  for width in range(limit):\n",
    "    for bits in helper.bitprod(width):\n",
    "      u = ops.Identity()\n",
    "      for bit in bits:\n",
    "        u = u @ base[bit]\n",
    "      gate_list.append(u)\n",
    "  return gate_list\n",
    "\n",
    "def create_unitaries_mitm(base, width):\n",
    "\n",
    "  gate_list = []\n",
    "  for bits in helper.bitprod(width):\n",
    "    u = ops.Identity()\n",
    "    for bit in bits:\n",
    "      u = u @ base[bit]\n",
    "    gate_list.append(u)\n",
    "  return gate_list\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1efa34d8",
   "metadata": {},
   "source": [
    "to_su2将任意 2x2 单位矩阵转化为 $SU(2)$ 矩阵\n",
    "\n",
    "F_dist定义了矩阵之间的距离度量\n",
    "\n",
    "generate_su2用于生成随机的$SU(2)$ 矩阵\n",
    "\n",
    "create_unitaries用于生成所有 0-limit 位的 bit 串，如：0, 1, 00, 01, 10, ... ，用来构造SK搜索空间。 create_unitaries_mitm适用于中间相遇策略。\n",
    "\n",
    "4.初始估计方法\n",
    "\n",
    "SK算法作为一种递推算法，需要以一个初始估计作为递推的起点，下面是两种初始估计方法：暴力搜索和中间相遇方法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b35ffd21",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bruce(gate_list, u):\n",
    "\n",
    "  min_dist, min_u = 10, ops.Identity()\n",
    "  for gate in gate_list:\n",
    "    tr_dist = F_dist(gate, u)\n",
    "    if tr_dist < min_dist:\n",
    "      min_dist, min_u = tr_dist, gate\n",
    "  return min_u\n",
    "\n",
    "def mitm(G,u):\n",
    "  \n",
    "  min_dist, min_u = 10, ops.Identity()\n",
    "  for i in range(len(G)-1):\n",
    "    for w in G[i+1]:\n",
    "      for v in G[i]:\n",
    "        a=to_su2(np.dot(v,w))\n",
    "        tr_dist = F_dist(a, u)\n",
    "        if tr_dist < min_dist:\n",
    "           min_dist, min_u = tr_dist, a\n",
    "      for v in G[i+1]:\n",
    "        a=to_su2(np.dot(v,w))\n",
    "        tr_dist = F_dist(a, u)\n",
    "        if tr_dist < min_dist:\n",
    "           min_dist, min_u = tr_dist, a\n",
    "  return min_u"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7e7ee19",
   "metadata": {},
   "source": [
    "暴力搜索直接线性搜索gate_list中的最接近项，而中间相遇（Meet-in-the-Middle）则是通过双向搜索逼近加速查找。\n",
    "\n",
    "5.GC分解\n",
    "\n",
    "SK算法的递归关键在于使用了一个所谓的GC分解的技术，通过两个操作的来回抵消，将误差压缩至更低阶。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ef1abc64",
   "metadata": {},
   "outputs": [],
   "source": [
    "def u_to_bloch(u):\n",
    "\n",
    "  angle = np.real(np.arccos((u[0, 0] + u[1, 1]) / 2))\n",
    "  sin = np.sin(angle)\n",
    "  if sin < 1e-10:\n",
    "    axis = [0, 0, 1]\n",
    "  else:\n",
    "    nx = (u[0, 1] + u[1, 0]) / (2j * sin)\n",
    "    ny = (u[0, 1] - u[1, 0]) / (2 * sin)\n",
    "    nz = (u[0, 0] - u[1, 1]) / (2j * sin)\n",
    "    axis = [nx, ny, nz]\n",
    "  return axis, 2 * angle\n",
    "\n",
    "\n",
    "def gc_decomp(u):\n",
    "\n",
    "  def diagonalize(u):\n",
    "    _, v = np.linalg.eig(u)\n",
    "    return ops.Operator(v)\n",
    "\n",
    "  axis, theta = u_to_bloch(u)\n",
    "\n",
    "  phi = 2.0 * np.arcsin(np.sqrt(np.sqrt((0.5 - 0.5 * np.cos(theta / 2)))))\n",
    "\n",
    "  v = ops.RotationX(phi)\n",
    "  if axis[2] > 0:\n",
    "    w = ops.RotationY(2 * np.pi - phi)\n",
    "  else:\n",
    "    w = ops.RotationY(phi)\n",
    "\n",
    "  ud = diagonalize(u)\n",
    "  vwvdwd = diagonalize(v @ w @ v.adjoint() @ w.adjoint())\n",
    "  s = ud @ vwvdwd.adjoint()\n",
    "\n",
    "  v_hat = s @ v @ s.adjoint()\n",
    "  w_hat = s @ w @ s.adjoint()\n",
    "  return v_hat, w_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87707dbc",
   "metadata": {},
   "source": [
    "6. Solovay-Kitaev 主递归算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "21662433",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sk_algo(u, gates, n, est):\n",
    "  \"\"\"Solovay-Kitaev Algorithm.\"\"\"\n",
    "\n",
    "  if n == 0:\n",
    "    if est == 'mitm':\n",
    "      return mitm(gates, u)\n",
    "    if est == 'bruce':\n",
    "      return bruce(gates, u)\n",
    "  else:\n",
    "    u_next = sk_algo(u, gates, n - 1, est)\n",
    "    v, w = gc_decomp(u @ u_next.adjoint())\n",
    "    v_next = sk_algo(v, gates, n - 1, est)\n",
    "    w_next = sk_algo(w, gates, n - 1, est)\n",
    "    return v_next @ w_next @ v_next.adjoint() @ w_next.adjoint() @ u_next"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65703b71",
   "metadata": {},
   "source": [
    "当 n=0 时，调用初始估计的搜索策略（bruce/mitm）\n",
    "\n",
    "否则：通过GC分解递归逼近误差项的 commutator，然后组合还原目标门。\n",
    "\n",
    "7.主程序"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3ee8cd02",
   "metadata": {},
   "outputs": [],
   "source": [
    "def main():\n",
    "  \n",
    "  num_experiments = 1           #随机生成SU2的个数\n",
    "  depth = 6                     #基本门序列长度\n",
    "  recursion = 3                 #SK算法深度\n",
    "  est='bruce'                   #初始估计方法\n",
    "\n",
    "  print('SK algorithm - est: {}, depth: {}, recursion: {}, experiments: {}'.\n",
    "        format(est, depth, recursion, num_experiments))\n",
    "  \n",
    "  base = [to_su2(ops.Hadamard()), to_su2(ops.Tgate())]\n",
    "\n",
    "  t1=time.time()\n",
    "\n",
    "  if est=='mitm':\n",
    "    depth = math.ceil(depth/2)+1\n",
    "    G = [create_unitaries_mitm(base, i) for i in range(depth)]\n",
    "  if est=='bruce':\n",
    "    G = create_unitaries(base, depth+1)\n",
    "\n",
    "  sum_dist = 0.0\n",
    "  for s in range(num_experiments):\n",
    "    u = generate_su2(s)\n",
    "    u_approx = sk_algo(u, G, recursion, est)\n",
    "    dist = F_dist(u, u_approx)\n",
    "    dist= math.log(dist,10)\n",
    "    sum_dist += dist\n",
    "    print('Distance: {:.6f}'. format(dist))\n",
    "\n",
    "  t2=time.time()\n",
    "  t=t2-t1\n",
    "  print('Mean time:{:.4f},  Mean Dist: {:.6f}'.\n",
    "        format(t/num_experiments, sum_dist / num_experiments))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66a23ae6",
   "metadata": {},
   "source": [
    "其中，num_experiments为随机生成SU2的个数，depth为基本门序列长度，recursion为SK算法深度，est为初始估计方法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "034f2daf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SK algorithm - est: bruce, depth: 6, recursion: 3, experiments: 1\n",
      "Distance: -1.464093\n",
      "Mean time:0.1329,  Mean Dist: -1.464093\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\97849\\AppData\\Local\\Temp\\ipykernel_17716\\3549546427.py:8: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return math.sqrt(1 - 0.5 * np.trace(u.adjoint() @ v))\n"
     ]
    }
   ],
   "source": [
    "if __name__ == '__main__':\n",
    "  main()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ea26476",
   "metadata": {},
   "source": [
    "运行main函数，打印出SK算法得到的门序列与目标门之间的距离（取log）以及算法运行的时间。\n",
    "\n",
    "## 3 参考文献\n",
    "[1] Dawson C M, Nielsen M A. The solovay-kitaev algorithm[J]. arXiv preprint quant-ph/0505030, 2005. https://arxiv.org/pdf/quant-ph/0505030\n",
    "\n",
    "[2] Amy M, Maslov D, Mosca M, et al. A meet-in-the-middle algorithm for fast synthesis of depth-optimal quantum circuits[J]. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2013, 32(6): 818-830. https://arxiv.org/pdf/1206.0758"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
